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There have been many discussions of two-mode models for Bose condensates in a double well 
potential, but few cases in which parameters for these models have been calculated for realistic 
situations. Recent experiments lead us to use the Gross-Pitaevskii equation to obtain optimum two- 
mode parameters. We find that by using the lowest symmetric and antisymmetric wavefunctions, it 
is possible to derive equations for a more exact two-mode model that provides for a variable tunneling 
rate depending on the instantaneous values of the number of atoms and phase differences. Especially 
for larger values of the nonlinear interaction term and larger barrier heights, results from this model 
produce better agreement with numerical solutions of the time-dependent Gross-Pitaevskii equation 
in ID and 3D, as compared with previous models with constant tunneling, and better agreement 
with experimental results for the tunneling oscillation frequency [Albiez et al., cond-mat/0411757]. 
We also show how this approach can be used to obtain modified equations for a second quantized 
version of the Bose double well problem. 

PACS numbers: 03.75.Lm,05.45.-a,03.75.Kk 
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I. INTRODUCTION 

The analogy between double Bose condensates, sep- 
arated by a barrier, and Josephson oscillations of su- 
perconductors [1] was apparently first suggested by Ja- 
vanainen [2] , and has been developed more thoroughly in 
a number of theoretical studies [3-18]. Especially from 
work in [4, 5, 10] and more recently in [12, 17], a rather 
elaborate picture of phase space dynamics has now been 
developed. The equations for tunneling oscillations of 
Bose condensates in a double well potential have been 
shown to resemble a pendulum whose length depends on 
the momentum. In the limit of small amplitude oscilla- 
tions, the equations are the same as for Josephson os- 
cillations for superconductors separated by a weak link 
[19]. It has also been shown that when atom-atom inter- 
actions exceed a critical value, the ensemble will remain 
trapped in one well while the phase continually increases, 
resembling a pendulum with sufficient energy to rotate. 

Experiments showing interference when condensates in 
a potential with a barrier were released [20] first stimu- 
lated interest in the problem of Bose condensates in a 
double well potential. More pertinent to the present dis- 
cussion are experiments that probe the evolution of the 
distribution between two or more wells of an optical lat- 
tice. Josephson oscillations have been observed in ID 
optical potential arrays [21]. Recently for a double well 
potential, both the regimes of tunneling and self-trapping 
of Rb atoms were observed [22]. In view of proposed 
extensions of these and other experimental techniques, 
[23-27], it seems appropriate now to reexamine the the- 
ory with the goal of developing models to deal with more 
diverse conditions. 

It is often assumed that the "tight-binding" approx- 
imation is valid, leading to what is known as the 
Bose-Hubbard model [28, 29], or discrete nonlinear 
Schrodinger equation [30, 31]. This model, which has 



been confirmed under the experimental conditions of [21], 
employs parameters for tunneling and on-site energy that 
are usually taken to be constant. One expects that with 
sufficiently large numbers of atoms, the atom-atom re- 
pulsion will cause the wavefunctions in a well to vary in 
size depending on the atom number, and consequently, 
the tunneling parameter and onsite energy might vary. 

We have found that it is possible to solve a more exact 
two-mode model based on symmetric and antisymmet- 
ric solutions of the Gross-Pitaevskii equation [11, 32]. 
For weak interactions, this new two-mode model pro- 
duces negligible differences from previous two-mode mod- 
els. However, for larger interactions, there are substan- 
tial differences and it turns out that the recent experi- 
ments [22] begin to sample the regime in which the dif- 
ferences are significant. In this report, we show that the 
new two-mode model implies a tunneling parameter that 
can vary with time, depending on the number and phase 
of the ensemble in each well, hence the name "variable 
tunneling model" (VTM). Despite the additional terms 
needed to produce this result, the equations eventually 
reduce to equations with the same form as the usual 
Bose- Josephson junction equations, but with parameters 
defined differently, and with one new term that can be 
significant in the case of strong interactions. Below, we 
compare results obtained with this model to those with 
a two-mode model with constant tunneling, with results 
of a multi-mode model, and then with numerical solu- 
tions of the time-dependent Gross-Pitaevskii equation 
(TDGPE). The parameters used in the two- or multi- 
mode models are obtained from numerical solution of the 
stationary Gross-Pitaevskii (GP) equation, so it is per- 
haps not surprising that the model that mimics the GP 
equation most closely also best reproduces results from 
the TDGPE. For very large interactions, results from any 
two- mode model will deviate from the TDGPE results, 
but agreement is most persistent with the VTM. Un- 
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der conditions of a particular experiment, effects of non- 
zero temperature and experimental uncertainties may be 
larger than the differences shown below. 

The present study is primarily limited to a mean-field 
approach using the GP equation, assuming that fluctua- 
tions and thermal excitations are negligible and without 
quantizing particle number. Because atom number is not 
quantized, the particle number difference and phase dif- 
ference of atoms in two wells or two modes are classical 
quantities in this approach. Considerable theoretical ef- 
fort has been devoted to the second quantized form [4, 7- 
9, f 7, 18]. We show below that our approach can be used 
to obtain more exact equations for quantization, and we 
apply these equations to the case of weakly interacting 
systems to show the connection between first and second 
quantized theories in this limited regime. As shown in 
[17], the classical patterns appear clearly in the quantum 
phase space picture with as few as 10 atoms. For experi- 
ments that involve on the order of 1,000 atoms, it seems 
useful to have an improved GP (mean-field) approach. 

An outline of this paper is as follows. Section 11 is 
devoted to ID models. We first (Section IIA) derive the 
new (VTM) two-mode equations, and then (Section IIB) 
compare with previous approaches. Section IIC lays out 
a multi-mode approach, Section IID discusses dynamics 
in phase space, Section HE gives equations for a second- 
quantized version. Experiments are of course in 3D, with 
some degree of transverse confinement. Therefore in Sec- 
tion III we present a formalism for 3D calculations and 
give a few results, including comparisons with the exper- 
imental results of [22]. 

II. MODELS 

A. Time-Dependent Gross-Pitaevskii Equation 

When the temperature is sufficiently low and when 
particle numbers are sufficient that second quantization 
effects are not important, the time-dependent Gross- 
Pitaevskii Equation (TDGPE) may be used for the wave 
function ip(x, t) for interacting Bose condensate atoms at 
zero temperature in an external potential V ex t{x). Let- 
ting h = m = 1, a dimensionless version is 

The relationship between g and g^o will be discussed 
in section IIIA. Here J dx\ip(x, t)\ 2 — N, where N is the 
number of atoms. Except in Sec. II F, A is not quantized, 
and the approach is strictly mean-field. 

We consider double- well potentials V ext (x) that are 
symmetric in x. Initially, we discuss one-dimensional 
versions. Under the above conditions, we will use results 
obtained with the TDGPE to test two-mode and multi- 
mode models discussed below. The TDGPE can tell us, 
for example, whether the phase is nearly constant as a 
function of x over an individual well. 



B. New Two Mode Model 

In many situations, a good approximation is obtained 
with a two mode representation of tp(x, t). In early work 
[4], wavefunctions localized in each well were used. Later, 
± combinations of symmetric and antisymmetric func- 
tions, as in [5, 10, 11] provided a more accurate formula- 
tion, and we follow the approach of [11] here: 

iP{x,t) = VN{Mt)<5>i(x)+Mt)<S>2(x)}; (2) 

•m.) = (3) 

where 

$±(x) = ±$±{-x); J dx^Qj = 5 itj ; i,j = +,-.(4) 

The $± will be assumed to be real, and to satisfy the 
stationary GP equations 

P±$± = -\ C ^ + V ext <$> ± +g\<5> ± \ 2 <5>± 1 (5) 

with g = gN. 

We can now define 

z(t) ee \Mt)\ 2 -\Mt)\ 2 , (6) 

4>(t) = 6 2 {t)-6 1 {t). (7) 

here 9i (t) are the phase arguments of the complex valued 
function tpi(t): ipi(t) = \tpi{t)\e l6 ^ tS> . The above normal- 
ization conventions lead to a constraint on the ipi(t): 

/oo 
^| 2 = A^|^(t)| 2 + |^W| 2 = l. (8) 
-oo 

Note that < &i($2) primarily occupies the left (right) well, 
but has nonzero density on the other side. In order to 
compare with results of the TDGPE, we define the the 
number of atoms in the left well as follows: 

r° n 

N L = dx\^(x,t)\ 2 = -+NzS, (9) 

J -oo A 

and we define S and An as 
S=[ dxQ+(x)\*-{x)\; An= — ~ Nr = 2zS.(10) 

J — oo 

From the ansatz (2), the TDGPE (1), and the GP 
equation (5), eventually one obtains differential equations 
for z and <j). We now briefly outline this derivation. In 
the following, these quantities will be used: 

Jij = 9 J $-(x)$ 2 (x) dx;(i,j =+,-) 

A 7 = 7__ - 7 ++ 
A/3 = /?_ - f3 + 



3 



A = 



107+- - 7++ - 7- 



B = (3- - (3+ + 7++ 2 7 " = A/3 - ^ 



C = 



7++ + 7— - 27+- 



Substitution of (2) and (5) into (1) yields 
.#a(t) 



z : ($+ + $ 



-($+ - $_ 



v 7 eft; 

= ^ [(^i(t) ± Mt)] [P± - gN\t> ± \ 2 } $ d 



5 7V 



2 ($ 3 ±p± + *|<f T g±) , (12) 



where 



P± = 2(^1 ± ^ 2 )-|'0i| 2 V'i T |-02| 2 V 2 ±^2+^1-01 
Q± = ±2(-0 2 -^ 1 ) + 5Vi|0i| 2 T 5^ 2 |V2| 2 ± V^-VlVv 

The usefulness of the <I>± basis is evident here, since in- 
tegrals with odd powers of <f> + or $_ vanish. From the 
above equations, including 4, the following equations for 
ipi,2(t) are obtained (fi, v = 1, 2; /1 7^ f): 

-^ + ^IVv| 2 + C^W, (13) 



In analogous coupled equations presented elsewhere, as 
in [4, 5, 10], the coefficient of ip v in the equation for i/> M 
is identified as the tunneling parameter, and it is usu- 
ally constant with time and independent of fi. In the 
above equation for ip^, there are additional terms in the 
coefficient of ip v that are functions of ipi{t) and ip2(t), 
hence varying with time. Since 1 1/^1 . 2 1 2 = (1± z )/2 and 
ipifa = (l/2)y / (l — z 2 )e lcl> , these extra terms depend on 
instantaneous values of both z and the phase difference, 
(p. For this reason, we will call this model the "variable 
tunneling model," or VTM. We will see that the addi- 
tional terms, although sometimes small, can bring this 
two-mode model into closer agreement with solutions of 
the TDGPE. 

Remarkably, despite some complexity of these addi- 
tional terms, relatively simple equations of familiar form 
can be obtained with no approximations beyond the as- 
sumption of a two mode representation of tp, as in (2). 
Equations (1), (2), (5), (12) and (3) are used. We obtain: 

dfi , Bz 

— = Az-\ = cos0 — Gzcos2</>; 



dt 

r j 7 

— = -By/l - z 2 sin 4> + C(l- z 2 ) sin 2(p. (14) 



These equations have the same form as analogous equa- 
tions in [10] except for the terms in C. They can be 
written in Hamiltonian form 



dH 



dH 
~dz~ 



(15) 



(11) with the Hamiltonian 



z 1 , 1 

H = A— — B^l-z 2 cos 4> + -C(l - z 2 ) cos 2<\> (16) 

This Hamiltonian is an integral of motion for a classi- 
cal system with generalized coordinates (z(t), (j>(t)) and 
dynamical properties (14) and will be referred later as 
a classical Hamiltonian. H is not equal to the expec- 
tation value ^jffl °f the quantum Hamiltonian H — 

— + V ext {x) + g\ip\ 2 within two mode approxima- 
tion (2). Since ip(x,t) defined as (2) is not an eigen- 

function of H, the expectation value ^^ff is not con- 
stant over time. However, the Hamiltonian (16) provides 
information about dynamics in phase space, including 
self-trapping, as will be discussed in section HE. 

In numerical work we often used a harmonic potential 
with Gaussian barrier of varying height and width: 



V ext {x) = \x 2 + V b e-^ 2 



(17) 



Equations (2) and (5) imply that distances are scaled 
by a x = \/h/Muj x , time by l/u x and energies, includ- 
ing Vb above, by hui x , where M is the atomic mass, and 
w x /2w is the harmonic frequency. This scaling will be 
used throughout this paper, and in particular in all the 
figures. To obtain numerical values for the overlap inte- 
grals -fij, where i,j = ±, we solved (5) using the DVR 
method [33, 34] with increasingly finer mesh, with iter- 
ations for each mesh to make the $± functions and the 
the nonlinear term self-consisent. Values for the 7^ are 
shown as a function of barrier height, Vb, for a=1.5, for 
gN=l, 10 and 100, in Fig. 1. For large enough Vb, 
all parameters 7^ are equal. As Vb decreases from the 
asymptotic region, 7++ decreases most rapidly because 
<f> + , with no node, is less "lumpy" than 

In Fig. 2, values for the parameters A, B, C and A7 
are shown. For small g, B is nearly constant and close to 
the value for a noninteracting gas, while A and C increase 
linearly with g. For higher values of Vb, larger values of 
g lead to more distortion of the C and A7 parameters. 
The parameter C is several orders of magnitude smaller 
than A, B and A7 except when g is large compared to 
one. When C is much smaller than A7, it is justified to 
neglect C but preserve the difference between B and A/3. 

If the term in C is neglected, we have effectively 
derived an alternative Bose-Josephson junction (BJJ) 
model, with revised parameters for tunneling phenom- 
ena between Bose condensates in a double well potential. 
It follows, for example, that the discussion of Joseph- 
son plasmons in [13] as Bogoliubov quasi-particles can 
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FIG. 1: Parameters 7++, 7 ,7-1 as a function of V& for gN 

= 1.0 (a); 10.0 (b); and 100.0 (c), calculated using (DVR) 
method applied to the stationary GP equation (5). The ver- 
tical arrows denote values of Vb for which /3+ = V&. 




-3 -2-1 01 
log, (gN) 



3-2-1 01 



FIG. 2: A, B, C, and A7 parameters as a function of g = gN, 
on a log-log scale. The three plots are for a = 1.5, and Vb 
= 4.0 (a); 7.0 (b); and 12.0 (c) in units of %uj z . The vertical 
arrows denote values of logg for which f} + — Vb- 



be applied to the present model, with appropriate sub- 
stitutions. However, when g is sufficiently large, the C 
term cannot be neglected. We emphasize that this term 
comes strictly from the nonlinear Gross-Pitaevskii equa- 
tion for Bose condensates in a double well potential and 
does not apply to superconducting Josephson junctions. 

A useful estimate for the condition for a two-mode 
model to be valid is that (3 + < With this rough 
criterion in mind, in many of the plots, we will denote 
this point by vertical arrows. By this test, in Fig. 1, two- 
mode models arc valid to the right of the arrows, while 
in Fig. 2, the regime of validity of two-mode models is 
to the left of the arrows. In reality, the transition is not 
sharp, as we will see below. 



C. Comparison with other model theories 

In [10], equations for ipi are derived with help of orthog- 
onal functions $1,2(2), defined as in (5) above. However, 
smaller terms were neglected: "This approximation cap- 
tures the dominant z dependence of the tunneling equa- 



tions coming from the scale factors ^1,2 ~ \/N\,2, but 
ignores shape changes in the wavcfunctions for N\(t) =/= 
iV 2 (f)." Effects due to shape changes were estimated to 
be small. We hnd this to be the case in certain regimes 
but not always. 

To make comparisons with the VTM, we write the 
equations from [10] taking % = m = 1 as above: 



.dtp 
<>—- 
dt 



- = (£? + C/i|^i| 2 W>i-^ 2 ; 



,; " 2 =(E% + UM 2 )ih~W 1 ; 



dt 



\Mt)\ 2 + \Mtt 



(18) 



where, for i = 1 or 2, 



dx 



Ui = gN J dx\ 



-\V^\ 2 + m 2 V ext 



I 1 



K 



I 



dx 



1 



(19) 



Since the tunneling term (K.) is constant with time, we 
will refer to this model as the "constant tunneling model" 
(CTM). In the CTM model, the functions $i, 2 are related 
to symmetric and antisymmetric functions <I>± as in (5). 
However the eigenvalues/chemical potentials (3± do not 
directly apply. To determine values for E® in terms of 
[3± and 7±,±, we introduce, for i = +, — , 



ti = j d.r 



1 d 2 ^ 
<P„ 

2 1 dx 2 



= &-.g(|$ 4 | 4 ) =A-7«- (20) 
For the quantities E° 2 introduced above, we obtain 
?0 _ jpO _ jp _ e + + e - 



E° 2 



E 



(21) 



Furthermore, 
U 1 = U 2 -- 



t/ = 5(|$i, 2 | 4 ) = ^|$+±$_| 4 ) 

6 7+ _+7__] = A + 2C. (22) 



= i [7++ 



In the symmetric/antisymmetric basis, the coupling term 
becomes 



K 



e+ _ A/3 - A7 



A7 
"4 



The Hamiltonian is 



Hctm = U Z -w- 2Wl - 2 2 



cos< 



(24) 



In part, the differences in the two approaches arise be- 
cause the wavefunctions, <I>i 2, extend somewhat into the 
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opposite well, as noted above. Rather than comparing 
the individual parameters, comparisons between the two 
two-mode models are better performed in terms of prop- 
erties that are independent of the model, and may be 
calculated also with the time-dependent Gross-Pitaevskii 
equation. One such property is the well-known Joseph- 
son plasma oscillation frequency [14], which is taken to 
be the oscillation frequency in the limit of small ampli- 
tudes of z and <f>. Another derived property is the onset 
of self-trapping at <f> — 0, which is usually labeled z c , the 
critical value of z. We will discuss this in subsection II E. 

In the limit of small z and <f), the equations for z and 
6 become 



CTM 



VTM : — 

dt 



dz 

~dt ~ 
= (2C- 



-2AC0; 



dt 



{U + 2K)z; 
= 2JC(U + 2JC): 



f t = { A + B-C)z; 
(B-2C)(A + B-C). (25) 



In every ID case we have considered, B — 2C > and 
A+B-C > 0. Numerical results obtained with the VTM 
and CTM models are shown in Fig. 3 in comparison 
with frequencies obtained with the TDGP equation. For 
j « 1, all three approaches agree well. For g=l and 
large Vb, the values for oj from the CTM are about 16% 
less than from the VTM, while for g=3, the asymptotic 
difference is about a factor of two. For larger values of 
g and for large \%, as illustrated for g=10 in Fig. 3c, 
JC becomes negative, hence u c becomes imaginary, and 
the real part of ujq plotted in Fig. 3 is zero. 

Values of B,2JC, and A7 for g = 10 and u z — 1.5 are 
shown in Fig. 4 (2K is plotted because (23) shows that 
in the limit 7^ — > 0, B — 2JC, and from (25), we see 
that B plays a role equivalent to 2JC). The region where 
JC < is clearly indicated. (3± are the actual eigenval- 
ues, which are calculated with the nonlinear interaction 
terms included. The quantities e± have no direct physi- 
cal meaning, so it is not surprising that they can lead to 
anomalous results. Note also that the putative regime of 
validity of two-mode models is to the right of the vertical 
arrows in each figure, and that for 5=10, JC is negative 
over most of this region. 

Thus from calculations of the Josephson plasma oscil- 
lation frequencies, we conclude that the additional terms 
derived in the VTM model take better account of non- 
linear interaction effects and produce better agreement 
with full TDGPE results. For low atomic numbers and 
weak interactions, these additional terms are not needed. 
It is also evident that as interactions increase in magni- 
tude, neither two-mode model reproduces TDGPE re- 
sults quantitatively. This will lead us to examine multi- 
mode models below. 

First, however, it will be helpful to take another per- 
spective by looking at results simply from the TDGPE. 
Fig. 5 shows |^(a;)| 2 and <p(x) as they evolve over one- 
half cycle under conditions in which (in a and b) the 
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FIG. 3: Comparisons of the oscillation frequencies for small 
z, (f> amplitude calculated from the CTM, the VTM and the 
TDGP equation, for gN = 1.0 (a); 3.0 (b); 10.0 (c). The 
insets in (a) and (b) show ratios of CTM and VTM results 
to TDGPE results. Vertical arrows denote values of Vb for 
which Vb — /?+. 



< 

0.2 



-B = p. - 13 + -Ay/2 
-2K = |3.-P + -Ay 
Ay = y_ - y ++ 



FIG. 4: Values for the parameters B,2JC, and A7 = 7 — 

7++ for g — 10 and values of the barrier height, Vb as indi- 
cated. Although JC becomes negative for Vb > 6, B remains 
positive. 



phase is nearly constant over each well, and (in c and 
d) with a larger g interaction parameter such that the 
phase over each well is not constant at a given time. In 
the latter case, the phase difference cannot be defined, 
and any two-mode model fails. 
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FIG. 5: (a) and (c) Evolution of ^(a;)! 2 over one-half cycle 
of tunneling oscillation, (b) and (d) Evolution of phase, <p(x) 
under the same conditions as in (a) and (c), respectively. The 
conditions for (a) and (b) are g = 2.0, Vb = 5.0; for (c) and 
(d), g = 10.0, Vb = 5.0, and in each case a = 1.5. In (a) and 
(c), the initial function ^(x)! 2 is denoted by thick solid lines. 
In (b) and (d), the initial value of the phase is everywhere 
zero. In (d) after the initial time, the phase is clearly not 
uniform over either well because of the strong interactions 
and low barrier. 



D. Multimode approximation 

From Fig. 3, we saw that there are deviations in the 
Josephson plasma frequency, ojq, between even the more 
exact (VTM) two-mode model and numerical solutions of 
the TDGPE. These deviations raise the question whether 
better agreement can be obtained by expanding the set 
of basis functions beyond simply <f> + and <&_. 

In this section we introduce a generalization of the 
VTM two-mode model. Starting from the TDGP equa- 
tion 



.dip 1 d 2 ip Tr , . , ,n , 



we introduce the following ansatz 



JV-l 



1p{ X ,t) = M*)« 



-i0 k t 



<Pk{x) 



(26) 



(27) 



k=0 



where <pk{x) satisfy the following equations 



$2s<p2s = - 
p2s+l<p2s+l 



2 dx 2 

1 d 2 4> 2s +i 

2 dx 2 



+ V ex t<p2s + g|0o| 2 02s 



+ Vext<p2s+1 
2, 



(28) 
(29) 



-g\<P\\ (P2 



a+l- 



Thus (pOA — 3>± as defined above, with normalization 
/ dxcpi{x)(pj{x) — Sij. Here we are effectively using the 
virtual excited states of the Gross-Pitacvskii equation 



rather than Bogoliubov quasi-particle states. Equilib- 
rium thermodynamics is not the goal here. Any orthonor- 
mal basis offers an extension of the two-mode model, and 
the quasi-particle basis is unnecessarily cumbersome for 
this application. Substituting the ansatz (27) into the 
GP equation (26) and using equations for <p2 S {x) and 
<p2s+i{ x ) with the orthogonality property, we obtain the 
following equation for the time depending amplitudes 
b r (t). 

ib r =-Eb2jJoo,2 3 ,re l( >-^ + ^)t 



-£&«+i7ii,2 J -+i,re < (- ft '+ 1+ ^>* 

3 

V b b*b -7 e i(-l3 n +fi m -l3s+l3 r ) 



(30) 



These are 2 J equation for real functions \bj(t)\ and 
arg(bj(t)), where J is number of modes. However there 
is the following constraint: ^2\bj(t)\ 2 = 1 which is a 

3 

consequence of the normalization condition for the wave 
function tp(x, t). Since also the overall phase is arbitrary, 
we effectively have 2 J — 2 equations for 2 J — 2 indepen- 
dent variables. Therefore, we define bj(t) = Cj(t)e laj ^\ 
with Cj(t) = \bj(t)\, and introduce the following variables 



A - r 2 - r 2 
(fir = Oi r —i — Ce r 



1, . . . , J — 1 
1,...,J-1 



It is not difficult to restate equations (30) in terms of the 
new variables. 

As in the case of VTM two-mode model, the main 
ingredients of multimode approximation are parameters 
Jklmn that can be found numerically from eigenfunctions 
of the Gross-Pitaevskii operator for the symmetric and 
antisymmetric "condensates." In making comparisons 
with two-mode model results and with numerical solu- 
tions of the TDGPE, we will use the number difference 
An(t) defined in (10), rather than z(t), which is not de- 
fined for the TDGPE. As an initial condition for the 
TDGPE, we use desired linear combinations of <I>± (re- 
labeled 0o,i in (28)). In a given experimental situation, 
the actual initial condition might differ and might need 
to be modeled more precisely. 

What our results show generally is that in circum- 
stances in which the VTM differs significantly from 
TDGPE, the time evolution curve is not sinusoidal, but 
is distorted by higher frequency components. Therefore 
one cannot easily extract a single frequency, for example, 
to correct the discrepancies exhibited in Fig. 3. Figure 6 
shows the actual time evolution curve for several cases. 
These curves should be viewed in light of the statement 
[31] that when /3+ < Vb, the tight-binding limit applies, 
or in our case, the two-mode VTM applies. As shown in 
this figure, the two- mode model agrees quite well with the 
TDGPE curve for g = 3.0, V b = 6.0, for which (3+ = 4.54 
is less than For larger g or smaller V b , the 2-mode 
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Time 



FIG. 6: Left column: The functions &+(x) and potentials 
V(x) for conditions indicated above each frame: g, Vt, with 
a z = 1.5 in each case. The last figure gives the chemical po- 
tential, /3+. Right column: time evolution of the fractional 
number difference, An times 10 3 for very small initial imbal- 
ance. The three curves are obtained from the 2-mode VTM 
model (long dashes), the 4-mode model (short dashes) and 
the TDGPE (solid curve). 



and TDGPE currves differ both in frequency and shape. 
In each of these cases, results obtained with a 4-mode 
model yield better agreement with the TDGPE curves. 
It is remarkable that this good agreement appears even 
for a very low barrier, Vf, = 2.0, for g=3.0. 



E. Phase space dynamics 

The evolution of z, <fi from the coupled equations, (14), 
closely resembles the dynamical evolution phenomena 
thoroughly discussed in [10]. We give a brief review to 
point out the differences arising from use of the VTM. 

To visualize the dynamics, it is helpful to view a plot 
of the Hamiltonian surface, H(z,<f>), as shown in Fig. 7 
for generic values of A, B and C. The surface is periodic 
in <f), with minima at z — 0, 4> = 2n7r and saddle points 
or maxima at z = 0, <p = (2n + l)7r, where n is an integer. 
Trajectories are horizontal curves (constant H) lying on 
this surface. 

Within cither two-mode model, self trapping occurs 




FIG. 7: Hamiltonian surface H(z,(j>) for Vo = 4, a = 1.5, 
g = 1. Trajectories lie on the surface, following contour lines 
that represent constant energy levels. 

for H above H s , the value of classical Hamiltonian at 
the saddle point. Critical values of z — z c are defined as 
values of z((j) = 2nir) = zq that give H(z, 4>) equal to H s . 
For \zq\ > z c , trajectories will not pass through z = and 
z will remain positive or negative. For the VTM model, 
the Hamiltonian given by (16), gives 

H s =H(0,tt) =B + ^ = H(z c ,0). (31) 

From this result and (16), we obtain 

z c y = J ^[B(A-B-C)] 1/2 . (32) 

For the CTM model, the Hamiltonian of (24) yields 
H s = H{p,ir) = 2JC = H(z c ,0), (33) 

so that 

z c , c = |[2/C(l/-2/C)] 1 / 2 . (34) 

Here the model breaks down when either K. < (see 
Fig. 4) or U - 2JC < 0. 

Before presenting results of calculations of z c , we recog- 
nize that as \z\ « | An| and g increase, as in Fig. 6, higher 
modes enter. The variation of An with time becomes ir- 
regular rather than close to sinusoidal, as shown by sev- 
eral plots obtained from calculations with the TDGPE 
in Fig. 8. For Fig. 8a and 8b (differing very slightly 
in An(0), but on opposite sides of An = z c ), closely re- 
semble results one would expect from a two-mode model. 
Figure 8c and d show irregular curves from the TDGPE 
in a regime where the two-mode model does not apply. 
In Fig. 8c, there are oscillations of An within the range 
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FIG. 8: Temporal evolution of An for various cases: (a) 
g = 3.0, V b = 9.0, z = An(0) = 0.0846; (b) g = 3.0, V b = 
9.0,An(0) = 0.0838; (c) g = 10.0, H = 5.5, An(0) = 0.411; 
(d) g = 3.0, V b = 4.0, An(0) = 0.735. 



An > before An eventually becomes < 0. Fig. 8d 
shows that An < is achieved for only a brief duration 
(between T = 25 and 29). Neither of these cases can be 
considered "self-trapping," but they are far removed from 
symmetric, periodic oscillations. Under such conditions 
of low barrier and/or strong interactions, it is somewhat 
arbitrary to make the distinction between self-trapping 
and not self-trapping. 

Nonetheless, we have attempted to establish criteria 
and apply them consistently so as to compare results from 
the CTM, VTM, and TDGPE approaches, as shown in 
Fig. 9. Here, z c values from (32) and (34) have been 
restated in terms of An c using (10) in order to compare 
with TDGPE results. For both g = 3.0 and 10.0, when 
Vb is high enough, there is good agreement between VTM 
and TDGPE results. CTM results are significantly lower 
for .g=3.0, while for #=10.0, as in Fig. 3, the fact that K 
becomes negative invalidates this approach in this regime 
of strong interactions. 

In the self-trapping regime, maximal and minimal val- 
ues of z(t) can be obtained by solving the equation z = 0: 



z = ± 



1 



1 - 



\B\ 



C-A 



± 



B 



C-A 



2H - A 

C-A 



.(35) 



Plus or minus signs in front of the square root correspond 
to different initial conditions for z (positive or negative 
respectively). An elegant discussion of dynamics, and 
separatrices, in phase space is given in [12] (explicitly for 
the case C=0). 

In [10], it was pointed out that closed trajectories on 
the surface of H can also occur around maxima on the 
lines (j) — (2n + l)n. These are the so-called 7r-phase 
modes. For the VTM, the condition for these maxima is 
that \B\ < \(A + C)\. The actual values, z„, at which 



An c (VTM):An c (TDGPE) 
V An c (CTM):An c (TDGPE) 




FIG. 9: Values for An c from the CTM, VTM and TDGPE 
models, for a) g = 3.0 and b) g = 10.0. The inset in (a) shows 
the ratio of CTM and VTM values to TDGPE results. 



these maxima occur can vary drastically from one model 
to the other. 

Even for the case of negligibly small overlap, the mo- 
mentum z(t) in the VTM differs drastically from the 
CTM when conditions place these two models on oppo- 
site sides of the transition to self-trapping. Far from the 
neighbourhood of the self trapping transition in g, a, Vn, 
the differences are less. 



F. Second Quantization 

Previous discussions of quantized versions of the Bosc 
double well problem [4, 9] were valid to first order in 
the overlap of the wavefunctions in each well. Also in 
the recent work [17], certain approximations are made 
for the wavefunction overlap. Now that we have a mech- 
anism for treating the overlap more exactly, there are 
new possibilities for extending the regime of validity of 
quantum approaches, which are necessarily based on two- 
mode models. 

The energy functional describing trapped BEC in 
terms of creation and annihilation operators tl>(x,t), 
ip^(x,t) can be written 



Ho — Hn 



Hn 



dx 



1 



^V 2 ^ + ^V ext tp 



(36) 



with the commutator [ip(x, t) 7 ^(x\ t)] = S(x — x'). 

As above, we will characterize the time evolution in 
terms of two modes that are predominantly (but not ex- 
clusively) located in the left and right wells. However, the 
derivation is easier when written in terms of the symmet- 
ric and antisymmetric functions, Q± rather than in terms 
of <j> 12 , because ($^_$ T ) = 0, whereas ($? 2 $2,i) 0- 
Wc therefore write a "mixed basis" expression: 



$ = -j= + + c 2 ($+ -*-)], 



in which 



ci, 2 = J dxtp(^ + ± 



(37) 



(38) 



are projections of ip. $± are solutions to the GP equation 
as above. In particular, the following form will be useful: 



2 



V ext * i =l3 i * i -gN\* i \ 2 \* i . 



(39) 



6ij. 



Also we have [c*, ct] 

Substituting 37 and 39 into the above equation for 77 , 
we obtain four terms: 



+ I c|c 2 + c\c-i 



\P+ ~ 7++ + /?- ~ 7— ] 
[/J+-7++-0-+7— ]} (40) 



Upon substituting 37 into the above equation for Hi, we 
obtain 16 terms, each with products of two creation and 
two annihilation operators, times integrals of the form 



(+)'( 



dx(<^+ +$_)*($+ -$_)7 (41) 



In particular 

Hi = c\c\c 1 c 1 (+) 4 + c\c\c 2 c 2 {-f (42) 
c\c\cic 2 + c\c\c 2 di + c\c\c\Ci + c\c\c\Ci (+) 3 (-) 



+ 



+ 



c\c\c\c 2 + c\c\c 2 ci + c\c\c\C\ 
+c\c\cic 2 + c\c\c 2 ci + c\c\c 2 c 2 



+ 



c\c\c\c 2 + c\c\c 2 c\ + c\c\c 2 c 2 + c\c\c 2 c 2 



(+) 2 (-) 2 
(+)(-) 3 



i=0 



Recalling the definitions 

All = 7 ++ + 67+- + 7 — 
A = U — 2C; A 7 = 7 __ 



4C* = 7 ++ +7__ 
-7++; A/3 = /3_ 



(43) 



27+-; 



we obtain 

(±) 4 = §/ dx($ + ±<i>_) 4 

1 2U 

= ^(7++ + 67+- +7— ) = ^; (45) 

(±) 3 ( T ) = § I -*-)(*l± 2*+*-+**) 

= ^(7 ++ -7-) = -^A 7 . 



(+) 2 H 2 = 9 -jdx(<S>\-&f 



1 2C 

= ^(7++ -2 7 +- +7— ) = ^p- (46) 

We wish to represent H 2 in terms of the following op- 
erators: 



N = N 1 +N 2 = c\c! + c\c 2 \ J x = -{c\c 2 - c\c x ) 



J y = ^(4 £ i - c\c 2 ); J z = -(4ci + c\c 2 ), (47) 



and Casimir element J 2 = ^ 



+ l), so that 
[Ji,Jj] = ie ijk J k . (48) 



Also wc will need 



A^ 2 2 

n 1 n 2 = — -j2. 



(49) 



Then the products of four annihilation and creation op- 
erators, the D(i) terms in (42), reduce to: 

A 2 



'Am + (4) 2 £i 

77(1) + 77(3) =4^-1)7,; 



77(0) + 77(4) = (c{) 2 c 2 + (clYci = — A + 2J 2 ; 



A 



77(2) = 4J 2 + 2AiA 2 = 4J 2 + — - 2J 2 .(50) 

Collecting terms, neglecting terms that are constant, we 
obtain 



77 2 = -J z A/3 + A7- 



2A 7 \ 4(A + C) 8C ? 

The quantum equations of motion read 

ji = *[H,Ji], (52) 

which yields 

Jx = (^0 + A 1 -^p-^J y -^-(J v J z + J z J y ) 
jy = -(a/3 + A 7 -^) Jx 

J^~(JxJz ~\~ J %Jx) ~ {JzJx JxJ z) 



4(A + C) 

•Jz — jy \' J y J z 1 < J z<Jy) 



(53) 



The above Hamiltonian, i/ 2 , is to be compared with 
expressions derived previously [4, 9, 12, 14, 17, 18]. 
Although a general second quantized Hamiltonian was 
written many years ago by nuclear physicists [35] (since 
known as the LMG model), most applications seem to 
involve simply the terms in J z and Jj. Using the op- 
erators defined above and assuming a symmetric double 
well potential, the expression in [14], for example, can be 
written: 



-£jj z + Ik j, 



(54) 



The comparison provides the following translation: 
2A 7 ^ ^ 8(A + C) 



£,, = A/3 + A 7 ■ 



N 



K 



N 



(55) 



The regimes defined in [14] then become (neglecting the 
2A 7 //V term: 



Rabi : 
Josephson : 
Fock : 



K 1 

--«-^i?«l; 



R 



8(A + C) 



1 K 

N << 8- J <<N 



A/3 + A 7 ' 
l«i?«]V 2 ; 



K 



N < R. 



(56) 



Thus for the second-quantized version as for the first- 
quantum GP equation version discussed above, we have 
obtained a Hamiltonian with a form similar to those pre- 
viously derived, but with slightly different parameters, 
and with extra terms that may be important for large 
atom-atom interactions. 

Our formulation provides a connection with experi- 
mental conditions through the Gross-Pitaevskii equation. 
The expectation value (2J X ) describes the difference be- 
tween the number of particles in the two modes, and is 
therefore an analog of the classical quantities, momen- 
tum, z(t), and number difference, An(t). The connec- 
tion between (2J X ) and An(t) can be most easily seen 
in the limit of very small interactions (small g), which 
is essentially the "Rabi regime" as defined in [14] and 
above. The following conclusions are based simply on an 
empirical evaluation of numerical results. 

For g < 10~ 2 , there are clear tunneling oscillations 
with frequency A/3 from the first term in H2 (A7 <C A/3 
here). These oscillations arc modulated by effects from 
the second term (in J 2 ) in H2, which increase with g. For 
constant g = gN, these modulations are independent of 
N over a large range of N, but undergo an additional 
modulation whose period decreases with N, as shown in 
Fig. 10. This suggests that there are various orders of 
time-dependent perturbations by which the second term 
in Hi perturbs the effect of the first. However, we have 
not been able to produce a quantitative perturbation- 
theoretic model. For long enough times, one observes 
the collapse and revival effects noted in [4] and shown in 
Fig. 11. For larger values of g, these structures no longer 
persist. 
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FIG. 10: In the limit of small g and small times, the temporal 
behavior of 2{J X ) and An(t) approach each other. All of these 
plots are for a = 1.5, 14=6.0. (a) g = 1.0 xl0~ 4 , R = 0.017; 
(b-d) g = 1.0 x 10~ 3 , R = 0.17. For (a) and (b), the quan- 
tum and classical periods are nearly the same, but (b) begins 
to show a decrease of oscillation amplitude in the quantum 
case, due to the second term in H-2- This modulation of the 
oscillations is shown over longer time in (c), for iV=100, and 
in (d), for iV=20. 



Since it has been difficult to experimentally observe 
even a single oscillation, and since it is difficult to con- 
trol N and the initial imbalance, these oscillatory pat- 
terns may be impossible to observe. We do find however, 
that the collapse and revival structure can persist even 
if there is a spread of initial values of N\ — Ni- Figure 
11 compares results for fixed iV=30, for an initial sharp 
Ni distribution, Ni — 27, with results for an initial flat 
distribution over the range 24 < Ni < 30. The stabil- 
ity of a part of the revival structure may occur because 
two parameters are needed to characterize a point on 
the sphere that is isomorphic to the SU (2) algebra used 
above. More extensive numerical results of phase space 
oscillations are given in [18], and detailed studies of av- 
erages in phase space using the Husimi distribution have 
been presented in [17]. 



11 



40000 




10000 _, £0000 

Time 



30000 



40O0O 



FIG. 11: Collapse and revivals for g=3.0 xl0~ 3 , i?=0.50, for 
precisely defined N = 30, but in (a) for precisely defined 
Ni = 27 and in (b) for a fiat distribution 24 < JVi < 30. 



III. CALCULATIONS IN 3D 

A. General Formalism 

In comparing with experimental results, the transverse 
confinement enters. In this study, we consider moder- 
ate transverse confinement, not approaching the Tonks- 
Girardeau regime [36]. We have extended the above 
methods to 3D as follows. We write the TDGPE first 
in MKS units, denoted by overbars: 



if). 



dt 



-V 2 + 



2m 2 ^ 

+V B + g 3D \$(5t,t)\ 2 ] ^(x,t) 



(57) 



where m is the atomic mass, and g 3 r> — 4irfr 2 a 3 D/m, a%D 
is the 3D scattering length, and j dx\tp(5i)\ 2 = N. The 
external potentials of interest here will include a purely 
harmonic term as given above, plus a barrier term as 
a function of z that will be chosen to be Gaussian or 
proportional to a cos 2 function, as in the experiments of 
[22]. 
We let 



y 



(58) 



and scale the coordinates and time as 

Xi = aiX t ; a 2 — fi/muJi] t = t/w z . (59) 
Then since 

J dx\yj(x)\ 2 =N = Jd5i\iP{5i)\ 2 = a x a v a z J dx|^(x)| 2 (60) 
%jj — (u x a.yU z )~ 1 l 2 '^) 1 and (57) becomes 



.&ip(x,y,z,t) 

1 — m — 



H{X,ip)ip(x,y,z,t) : 



W(A,V) 



+ y +V B +4w V (^j mx,y,z,t)\ 2 (61) 



where A represents the arguments r), Vb, N, a^o and a z . 
An ansatz analogous to (2) can now be introduced: 



V>(x,i) 



JV[V>i(*)*i(x)+V2(*)*2(x)]; (62) 



$ l,2(x) 



$+(x) ± $_(x) 

71 ' 



(63) 



where 



®±(x,y,z) =±$±(x,y,-z); / dxdydz$% = 1. (64) 



(65) 



The stationary GP equations take the form: 
0±®±{x,y,z) = H(\,®±)$±(x,y,z). 



Because the transverse wavefunction is very nearly 
Gaussian, some authors have simply assumed a Gaus- 
sian, possibly with a z-dependent width, and obtained 
modified equations [37, 38] for what we have called tpi(z). 
Because we wanted to consider cases where the Gaussian 
form may not be valid, we used general 3D algorithms. 
Initial $± wavefunctions were obtained by diagonalizing 
the DVR Hamiltonian [34] using sparse matrix techniqes 
[39], which made calculations with >120,000 mesh points 
possible in minutes on a PC. To calculate the time evo- 
lution, the split-operator method [40] with Fast Fourier 
transform [41] was used, requiring an hour or more of f«2 
GHz CPU time, in view of the small time steps required. 

From the <3?± functions calculated from the 3D time- 
independent Gross-Pitaevskii equation, one can also ob- 
tain the parameters /3±,7i,j, A, B and C as in Section II, 
to provide a two-mode representation of tunneling oscil- 
lations in 3D. In translating results from ID to 3D for 
9id = 93D = ^V a 3D/ct z , we find that, in the limit of 
weak interactions and r\ > 1, the 7^ functions for 3D are 
a factor of 2-7T smaller than the corresponding ID 7^ func- 
tions. The explanation touches on the basic properties 
of tight transverse confinement. 

If the transverse confinement is symmetric in x and 
y and is tight enough, the 3D wavefunction $ + (x,y, z) 
can be factored into a function of z and a function of 
P 



\J (x 1 + y 2 ). Then if also the interactions are weak 



enough, the p function will be a Gaussian: 



$+{x,y,z) 



<M P ,z) = e -" 2 / 2 <Mz). 



(66) 



The normalization condition is 



J dxdydz$ + (x,y, z)' 
= 2n J dppe- p2 J dzM z ) 2 = k j dzM z ) 2 - (67) 
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FIG. 12: Values for 

lijihj — +> — Irom 3D calculations with 
77=1 and 100 (left axis), as compared with a ID calculation 
(right axis). Solid lines denote 7++, short dashes 7h — , and 

longer dashes 7 For 3D, the 7^ are approximately 2n 

larger than for ID. 



Under the above conditions, and if g\ d — g3D, then to 
within a constant of proportionality, 4>\{z) will also be 
a solution of the ID problem: <fri(z) oc $^(2). For the 
ID problem, J dz\^ 1 ^ D {z)\ 2 — 1, so from the different 
normalizations, we see that, under all the above stated 
conditions, 



(68) 



Mz) = ^] D (z). 

The 3D version of j ++ becomes 

7++ = / dxdydz\$ 3 + D \ 4 

= 2nJ pdpe- 2 ^ J dz^izf = 2n±?±± = 2±±. (69) 

Similar relations hold also for 7^ and 7+^. For larger 
interactions, the p dependence is not exactly Gaussian, 
the functions Q^P no longer factorize, and the parameters 
7?° deviate from the above relations. Figure 12 shows 

plots of 7 ++ ,7__ and 7^ from 3D calculations with 77 

= 1 and 100, as compared with ID results. For 77=100, 
the wavefunction is more concentrated than for 77 = 1, so 
the values for 7^ are slightly larger. Each is 5 to 6 times 
smaller than for the 3D case. Otherwise the dependences 
on Vb are very similar. 

There are other differences between 3D and ID prop- 
erties. The difference energy, A/3, and hence also the pa- 
rameter B decrease more rapidly as a function of barrier 
height. Figure 13 compares the parameters A, B and C 
in 3D (77 = 1) and in ID, for the case g — 10. Evidently, 
finite transverse confinement decreases the difference be- 
tween the symmetric and antisymmetric condensate en- 
ergies. The differences are much the same for 77 =100 as 
for 77=1. Also for 77 = 1 g—10, Fig. 13b shows that the 
plasma oscillation frequency in the limit of small z, <j>, for 
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FIG. 13: (a) A comparison of calculated values for the pa- 
rameters A, B and C for ID and for 3D, for r\ = 1 and <?=10. 
(b) A comparison of Josephson plasma oscillation frequencies, 
u>o, for ID and 3D, also for 77 = 1 and <7=10. 



barrier height, VJ, > 5, is even less than a factor of 2ir 
smaller in 3D than in ID. This statement has been found 
to be true for 77=1 and 100, and g up to 10. 

We conclude that two-mode models tend to be even 
more valid in 3D than in ID. 

Linear combinations of <i>± functions provide the initial 
condition for the TDGPE, for which we use the split op- 
erator approach with fast Fourier transform [40, 41]. To 
be able to compare TDGPE results with (25), we use a 
very small initial imbalance (z =0.002) for the TDGPE 
calculations. For the two-mode models, parameters are 
obtained from wavefunctions calculated with the time- 
independent 3D GP equation, as for ID results above. 
The results for ui are shown in Fig. 14a-c. (Fig. 14d per- 
tains to the experiments of [22] as discussed below). The 
plasma oscillation frequency obtained from the TDGPE 
increases rapidly beyond g — gN w 3. The two-mode 
model results match the TDGPE results well for g < 1. 
Fig. 14a, for Gaussian barrier of height Vb = 5Tllu z , 
shows good agreement for both the VTM and CTM with 
TDGPE results, up to g = 100. On the other hand, 
when the barrier height is raised to 8huj z , the CTM fails 
for g > 30, for both 77 = 1 (b) and 77 = 100 (c). For 
the latter, the VTM result begins to deviate significantly 
from the TDGPE value around 5=100. The failure of the 
CTM here is analogous to the situation shown in Fig. 3, 
and occurs because JC becomes negative. 
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FIG. 14: Josephson plasma oscillation frequencies calculated 
from the TDGPE and from VTM and CTM models, using 
(25) and the /3±,e± and y± parameters obtained from the 3D 
wavef unctions, $±. For (a-c), Gaussian barriers were used, as 
all the ID calculations. For (d), a cos 2 function to reproduce 
the barrier in [22]. Other parameters were (a) r/=1.0; 14=5.0; 
(b) t/ = 1.0, V b = 8.0; (c) r\ = 100, V b = 8.0; (d) Correspond- 
ing closely to the experiments in [22] (see text), v x — 66 Hz, 
v y = 90 Hz, v z = 78 Hz, barrier Vb = 5.28 Tiw z times the 
cos 2 function given in the text. The values plotted are for 
«o = 0.28 as in the experiments, rather than for the limiting 
case of small z,<f>, hence are labeled u> rather than loq. 



B. Comparison with Recent Experiments 

Very recently, the first quantitative experimental ob- 
servations of oscillations of Bosc condensates in a dou- 
ble well potential have been performed [22]. The har- 
monic confinement was created by overlapping tightly 
focussed Gaussian laser beams. The harmonic frequen- 
cies were 66, 90 and 78 Hz in what we will call the x,y 
and z directions. To produce the double well, an op- 
tical standing wave from two beams of wavelength 811 
nm, at an angle of 9° were added, producing a barrier 
of the form Vj, cos(ztt / w) 2 , with VJ, = 413(20) Hz, and w 
= 5.20(20)/im. 1150 87 Rb atoms were loaded into this 
trap. We have modeled this experimental configuration 
and find the effective value of gN — 58.8 from 61, which 
corresponds to g ss 10 in ID simulations. 

The reported experimental period of oscillation for 
z = 0.28(6) was 40(2) msec. [22], which corresponds to 
the value indicated by the large diamond in Fig. (14)d. 
We obtain values very close to this observed tunneling 
frequency with both the VTM and TDGPE by using a 
value for a^D = 100ao (where ao is the Bohr radius) [42]. 
For this initial value zq, (although not for very small 
values of zq), self-trapping occurs with the CTM model, 
when based on solutions of the Gross-Pitaevskii equation. 
The calculated CTM frequency drops rapidly before this 
point, as shown in Fig. (14)d. 



Note that by reference to Figs. 14a-c, we conclude that 
as long as the temperature is sufficiently low (the temper- 
ature for the experiments of [22] was immeasurably low), 
the aspect ratio is not important, as results for 77 = 100 
and =1 are very similar, but the relatively large value of 
the nonlinear term is important in determing the validity 
of two-mode models. 

Using the TDGPE, we have calculated a value of 
z c =0.39 for the stated conditions of these experiments, 
which is consistent with the observed oscillations at zq — 
0.28(6) and self-trapping for z(0) = 0.62(6), but lower 
than the value of z c = 0.50(5) quoted in [22]. In this pa- 
per, the authors performed calculations with the trans- 
verse Gaussian model of [37] and obtained good agree- 
ment with experimental observations. Our contribution 
is simply to show that a two-mode model, with param- 
eters from GP eigenfunctions, also comes quite close to 
reproducing the experimental observations. 

The other experiments that helped to motivate this 
study were performed at NIST, MD, with 87 Rb atoms 
in a "pattern-loaded" optical lattice. The atoms were 
first loaded into a coarse lattice from Bragg-diffracted 
laser beams, and then a finer lattice was turned on, 
such that every third lattice site was occupied [23]. We 
are presently working to develop a modification of the 
present approach to deal with such phenomena in a pe- 
riodic lattice. 



IV. CONCLUSIONS 

By rigorous solution of coupled equations for the sym- 
metric and anti-symmetric wavefunctions for a Bose con- 
densate in a double well potential, we have derived equa- 
tions for a new two-mode model that provides for vari- 
ation of the tunneling parameter with time, depend- 
ing on the differences of number and phase of atoms 
in the two wells. We have compared results from this 
"variable tunneling model" (VTM) with results from 
other two-mode models, from multi-mode models that 
we have constructed, and from solutions of the time- 
dependent Gross-Pitaevskii equation (TDGPE). In mak- 
ing these comparisons, we numerically compute wave- 
functions from the stationary Gross-Pitaevskii equation 
and use appropriate integrals of these wavefunctions in 
the model equations. For small values of the non-linear 
interaction term and moderate potential barriers, all 
the models agree nicely. When the nonlinear interac- 
tion term exceeds a certain value, the tunneling parame- 
ter in the usual "constant tunneling model" (CTM) be- 
comes negative, and thus the Josephson plasma oscilla- 
tion frequency becomes imaginary. The VTM remedies 
this problem, and produces better agreement with re- 
sults of the TDGPE. We have performed such compar- 
isons for ID situations and also for 3D situations, for 
which we have obtained 3D solutions of the stationary 
and time-dependent GP equations. The recent experi- 
mental observations of tunneling oscillations and macro- 
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scopic self-trapping of [22] are in the regime of moder- 
ately strong non-linear interactions because of the large 
number of atoms (1150 87 Rb atoms). Results from both 
the TDGPE and VTM for the observed tunnelling os- 
cillation frequency are in good agreement with the ex- 
perimental value. However, under the conditions of the 
experiment, the CTM when based on parameters from 
the GP equation, leads to self-trapping rather than oscil- 
lation. 

We also have applied our approach to derive an im- 
proved Hamiltonian for quantum calculations, but find 
no reliable standards to compare this approach with 



other quantization approaches. What we have not inves- 
tigated here are damping effects of thermal excitations, 
as considered in [6] and [43]. 
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